Discovery of a Dysprosium Metallocene Single-Molecule Magnet with Two High-Temperature Orbach Processes

Magnetic bistability in single-molecule magnets (SMMs) is a potential basis for new types of nanoscale information storage material. The standard model for thermally activated relaxation of the magnetization in SMMs is based on the occurrence of a single Orbach process. Here, we show that incorporating a phosphorus atom into the framework of the dysprosium metallocene [(CpiPr5)Dy(CpPEt4)]+[B(C6F5)4]− (CpiPr5 is penta-isopropylcyclopentadienyl, CpPEt4 is tetraethylphospholyl) leads to the occurrence of two distinct high-temperature Orbach processes, with energy barriers of 1410(10) and 747(7) cm–1, respectively. These barriers provide experimental evidence for two different spin–phonon coupling regimes, which we explain with the aid of ab initio calculations. The strong and highly axial crystal field in this SMM also allows magnetic hysteresis to be observed up to 70 K, using a scan rate of 25 Oe s–1. In characterizing this SMM, we show that a conventional Debye model and consideration of rotational contributions to the spin–phonon interaction are insufficient to explain the observed phenomena.

Alternative synthesis leading to a mixture of [2][B(C6F5)4] and [3] in cold hexane (-40 °C, 10mL) with a glass-coated stirrer bar. The reaction was warmed to room temperature. The resulting suspension was sonicated for 50 minutes and then stirred for 48 hours, which produced a yellow powder. After letting the solid settle, the hexane was decanted away and the residue washed with hexane (5 × 10 mL). The resulting yellow powder was dried under vacuum. Toluene (5 ml) was added and the suspension was heated to 90°C to give a yellow solution with an orange oil. Cooling to room temperature and storing for two days produced bright yellow crystals. The solvent was decanted away and the crystals were washed with cold toluene (3 × 5 mL) and hexane (3 × 5 mL), giving a mixture of crystals of [2][B(C6F5)4] and [3][B(C6F5)4] as the major and minor product, respectively. S9 X-ray crystallography X-ray diffraction data for [(Cp iPr5 )Dy(Cp Et4P )(BH4)] (1) and [3][B(C6F5)4] were collected on an Agilent Gemini Ultra diffractometer at the University of Sussex. Data for [2][B(C6F5)4] were collected using a Rigaku microfocus rotating anode diffractometer at the University of Sussex using Cu-K radiation. Data for [2][B(C6F5)4] at 100 K and 30 K were also collected on a Bruker Smart APEX II diffractometer, using Mo-K radiation. Structures were solved in Olex2 with SHELXT using intrinsic phasing and were refined with SHELXL using least squares minimisation. [4][5][6] Anisotropic thermal parameters were used for the non-hydrogen atoms and isotropic parameters for the hydrogen atoms. Hydrogen atoms on carbons were added geometrically and refined using a riding model. Hydrogen atoms for BH4 fragments in compound 1 were located from the Q peaks around the boron atoms and isotropically refined. Structures were deposited at the Cambridge Structural Database and assigned the references codes shown in Table S1. S10 Table S1. Crystal data and structure refinement for 1, [2][B(C6F5)4] (on two instruments at 100 K, and at 30 K on one) and [3][B(C6F5)4].  (2) Cpc-Dy-Cpc 147.69(4) S14  Table S3 and is different to the atom labelling scheme in the corresponding CIF.  Figure S7. Molecular structure of the cation 2. Collected on a Bruker Smart APEX II diffractometer at 100 K. Thermal ellipsoid representation (50% probability). Green = dysprosium, orange = phosphorus, black = carbon. Hydrogen atoms are omitted for clarity. S17 Figure S8. Molecular structure of the cation 2. Collected on a Bruker Smart APEX II diffractometer at 30 K. Thermal ellipsoid representation (50% probability). Green = dysprosium, orange = phosphorus, black = carbon. Hydrogen atoms are omitted for clarity. Figure S9. Thermal ellipsoid representation (50% probability) of the structure of the cation 3. Green = dysprosium, orange = phosphorus, black = carbon, pink = boron, white = hydrogen (only non-C-H hydrogen atoms are shown).
Powder X-ray Diffraction Powder XRD data were collected using synchrotron radiation at the Diamond Light Source, beamline I11, at room temperature and 100 K. All experiments were conducted by Dr Sarah Day. The sample took the form of crushed polycrystalline material and was measured in a glass capillary. A beam current of 300 mA was used and data files were obtained using the MAC detection system. Instrument calibrations were performed using Si powder standard (NIST SRM640c). The wavelength was  = 0.826838(5) Å and the 2 zero-point was 0.00159°. Magnetic property measurements Samples were prepared by adding the crushed crystalline materials and eicosane into 7mm NMR tubes. The tubes were flame sealed under a static vacuum. The eicosane was melted in a water bath at 40°C to prevent crystallite torquing. The direct current (D.C.) magnetic susceptibility and magnetization data were collected using a Quantum Design MPMS3-VSM SQUID magnetometer in cooling mode. Alternating current (A.C.) magnetic susceptibility measurements were performed using a Quantum Design MPMS XL-7 SQUID magnetometer using an oscillating field of 5 Oe. Measurements of the magnetization decay were performed by first magnetizing the sample in a field of 7 T, and then removing the field and measuring the magnetization over time with the VSM option. Diamagnetic corrections were performed using Pascal's coefficients.    Figure S19. Temperature dependence of the isothermal and adiabatic susceptibility obtained from fitting the frequency-dependent susceptibility data (equation S1 and S2) using a single Cole-Cole model. The temperature dependence suggests two processes with a relative magnitude of 35% of the two relaxation processes. Figure S20. Frequency-dependent AC susceptibility data taken at 75 K, clearly demonstrating two relaxation processes. When fitted to two independent Cole-Cole processes (dashed lines), the resultant fit (solid lines) gives a ratio of the isothermal susceptibility of 38% (0.137 vs 0.052 cm 3 mol -1 ), strongly suggesting two processes. The data therefore cannot be fitted with two interdependent Cole-Cole processes. Figure S21. Cole-Cole plots for the AC susceptibilities in zero DC field for [2][B(C6F5)4] from 105-119 K. Solid lines represent fits to the data using equations 1 and 2, which describe ' and '' in terms of frequency, isothermal susceptibility (T), adiabatic susceptibility (S), relaxation time (), and a variable representing the distribution of relaxation times ().

S72
Computational details Two different types of geometry optimizations were carried out: a full optimization of the cation 2 and an optimization of the positions of the hydrogen atoms with the coordinates of heavier atoms frozen to their crystal-structure positions of the cations 2 and 3. Both the major and minor disordered components of 2 and 3 were considered in the calculations. We will refer to the structures of the major components with the optimized hydrogen positions as 2a and 3a, and to the minor components as 2b and 3b when distinction between the structures is computationally relevant. The geometry optimizations were carried out using density functional theory (DFT).
The DFT calculations were conducted using the ADF 2019 software [1]. Scalar relativistic effects were accounted for with the zeroth order regular approximation (ZORA). [2] The standard Slatertype basis sets designed for ZORA calculations were used throughout. [3] A valence triple-ζ quality basis set with two sets of polarization functions (TZ2P) was used for the Dy 3+ ion, valence triple-ζ quality basis set with a single set of polarization functions (TZP) was used for the P atoms and the C atoms in the Cp rings, and a valence double-ζ quality basis with a single set of polarization functions was used for the remaining C atoms and the H atoms.
The geometry optimization were carried out using the pure PBE exchange-correlation (XC) functional [4] and Grimme's DFT-D3 empirical dispersion correction [5] in conjunction with the Becke-Johnson (BS) damping function [6]. The "NumericalQuality" keyword in ADF was set to "Good" to avoid convergence issues due to numerical noise. To simulate the static correlation effects within the open 4f shell, the two unpaired β electrons were equally distributed over the seven 4f orbitals to yield non-integer occupation numbers. The geometry convergence tolerances were increased to 10 -4 , 10 -4 , 10 -3 and 10 -1 atomic units for the energy, gradient, bond length step size and bond angle step size, respectively. In case of the full optimization, a further frequency calculation was carried out to ensure that the stationary geometry corresponded to a minimum on the potential energy surface and to provide the normal modes needed for the spin-phonon calculations (vide infra). To smoothen convergence and to simulate static electron correlation effects, the nine 4f electrons were equally distributed over the seven 4f orbitals. In practice this means, that the occupations of the seven highest-energy β orbitals are fixed to the fractional values 2/7 ≈ 0.285714.
The hyperfine coupling tensors were calculated using geometry 2a. The calculation was carried out using the PBE0 hybrid XC functional. [4,7] The basis sets of the target atoms (P and Dy) were augmented to TZ2P-J. [8] The TZ2P-J basis corresponds to the TZ2P basis but has additional steep s functions to better account for the electron density near the nucleus. The "NumericalAccuracy" keyword in ADF was increased to "VeryGood". The hyperfine coupling tensors were calculated using the ZORA methodology implemented in ADF with perturbative treatment of spin-orbit coupling (SOC) effects. [9] The magnetic properties were calculated using multireference electron correlation methods as implemented in the Molcas quantum chemistry software version 8.4 [10] and OpenMolcas version 21.02 [11]. The static properties on structures 2a, 2b, 3a and 3b were calculated using OpenMolcas whereas the spin-phonon coupling constants were calculated with Molcas. First, a set of stateaveraged (SA) complete active space self-consistent field (CASSCF) calculations [12] were carried out. The active space consisted of the nine 4f electrons in the seven 4f orbitals. All 21 sextet, 224 quartet and 490 doublet states were solved in three separate state-averaged calculations. Spinorbit coupling (SOC) was accounted for using the restricted active space state-interaction (SO-RASSI) approach [13] where the SOC operator was constructed in a basis of the CASSCF eigenstates. All SA-CASSCF states up to an energy-cutoff of 50,000 cm -1 were included in the SO-RASSI treatment; this set of states consisted of all 21 sextet, lowest 128 quartet and lowest 130 doublet states. The static magnetic properties, the ab initio CF decomposition and the effective relaxation barrier were calculated using the SINGLE_ANISO routine in Molcas and OpenMolcas. [14,15] Scalar relativistic effects were included using the scalar version of the exact two-component (X2C) transformation as implemented in Molcas. [16] The atomic mean-field integral (AMFI) formalism was used in the construction of the SOC operator. [17] Relativistic atomic natural orbital (ANO-RCC) basis sets were used in all multireference calculations. [18] A polarized S73 valence quadruple-ζ basis set (VQZP) was used for the Dy ion, polarized valence triple-ζ basis sets (VTZP) were used for the P atom and the C atoms in the Cp rings, whereas a polarized double-ζ basis sets (VDZP) were used for the remaining C atoms and a simple double-ζ basis (VDZ) was used for the remaining H atoms. Cholesky decomposition with a threshold of 10 -8 atomic units was used to reduce the necessary storage space for the two-electron integrals.
The spin-phonon coupling constants of phonon α between CF eigenstates and were calculated as the matrix elements of the derivative of the CF operator ^ with respect to unitless normal modes : From the Hermicity of ^ it follows that ∨ ∨. The CF operator is written in terms of the equivalent operators ^ ^ of rank k and component q [19] where is a complex CF parameter. The rank k has even, non-negative values. The largest value of k is 2J. The components range in integer steps from q = -k to k. The equivalent operators are used in the Iwahara-Chibotaru notation [20], where the matrix element between angular momentum eigenstates is given by a ratio of two Clebsch-Gordan coefficients: .
The phonons were approximated as the gas-phase normal modes of the molecule. This approximation neglects both phonon dispersion as well as all acoustic modes. The spin-phonon coupling constants were constructed from Cartesian derivatives of the CF operator by [21] = ℏ , ,

^ ^ ,
where the index i runs over Cartesian displacements of all atoms, are the normalized eigenvectors of the Hessian in mass-weighted coordinates, is the angular frequency of normal mode α, is the atomic mass of the atom displaced in displacement , the indices M and M' run over the projections of the angular momentum J, are the CF eigenvectors.
The normal modes were determined from the full DFT optimization of 2. The Cartesian derivatives of the CF parameters were calculated by conducting a series of multireference calculations on geometries displaced along each Cartesian displacement direction. Each atom was displaced in five steps of 0.001 Å in positive and negative directions along each three Cartesian coordinates. Thus, a totality of 2671 displacements (including the equilibrium structure) were calculated. The large number of calculations necessitated a reduction in the level of theory as compared to the calculation of the static magnetic properties. Only the 21 sextet states were considered in the SA-CASSCF and SO-RASSI calculations and electron correlation effects outside the 4f orbital space were neglected. A VTZP basis was used for the Dy 3+ ion, VDZP basis sets were used for the P atom and the C atoms in the Cp rings and a valence double-ζ basis (VDZ) sets were used for the remaining C atoms and the H atoms. The optimization of the equilibrium geometry introduces some deviation between the static magnetic properties calculated for the crystal-structure geometry and for the optimized equilibrium geometry, although all considered properties are qualitatively similar. To keep the main text consistent, the energies of the KDs are discussed in terms of the values calculated for the crystal structure geometry.

S74
The CF parameters were calculated for each displaced geometry using the ab initio CF theory. [15] The Cartesian derivatives of the parameters were extracted by conducting a linear fit on each parameter using the ten displaced geometries and the equilibrium geometry. The complex CF parameters are only determined up to a phase, which varies between the calculations on different displaced geometries. To eliminate the phase-dependence of the results, the phase of each CF parameter was fixed to that in the equilibrium geometry and only the magnitudes of the parameters were varied in the fit. The effect of this approximation on the quality of the results was estimated by diagonalizing the CF operator in each displaced geometry with and without the phase-shift and comparing the results. In each 2671 geometries both the maximum and mean error in each eigenvalue was less than 1 cm -1 . Thus, the effect of the phase-approximation on the results can be considered negligible. The numerical stability of the fits was studied by examining the coefficient of determination of each fit. In the case of large derivatives, the coefficient was usually larger than 0.98. However, in the case of small derivatives very low values of some coefficients was observed.
To eliminate derivatives which could not be distinguished from numerical noise, all Cartesian derivatives with absolute values smaller than 1 cm -1 Å -1 were removed from the results. This necessarily introduces some error to the quantitative values of the spin-phonon coupling constants. However, the values of coupling constants with absolute values considerably smaller than 1 cm -1 still largely depend on the Cartesian derivatives with values larger than 1 cm -1 Å -1 ; thus, the qualitative picture remains the same.